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Abstract — We present an Alternating Direction Method of 
Multipliers (ADMM) algoritlini for solving optimization prob- 
lems with an £i regularized least-squares cost function subject 
to recursive equality constraints. The considered optimization 
problem has applications in control, for example in £i reg- 
ularized MPC. The ADMM algorithm is easy to implement, 
converges fast to a solution of moderate accuracy, and enables 
separation of the optimization problem into sub-problems that 
may be solved in parallel. We show that the most costly step 
of the proposed ADMM algorithm is equivalent to solving an 
LQ regulator problem with an extra linear term in the cost 
function, a problem that can be solved efficiently using a Riccati 
recursion. We apply the ADMM algorithm to an example of 
£i regularized MPC. The numerical examples confirm fast 
convergence to moderate accuracy and a linear complexity in 
the MPC prediction horizon. 

I. Introduction 

In this paper we consider optimization problems with an 
£i regularized least-squares cost function subject to recursive 
equality constraints. This has applications in control. The 
least squares part is standard in this context and penalizes 
deviations of the states from the set-point at the same time 
as keeping the control signal small. The ^i-norm regular- 
ization of the cost function promotes sparse solutions, i.e. 
a solution with many zero entries, [1]. The cost function is 
known as LASSO, [2]. LASSO is a well-known method in 
statistics and machine learning, and it has gained a lot of 
interest in other research communities as well, e.g. system 
identification, [3]. 

We propose to solve the optimization problem using an 
algorithm called Alternating Direction Method of Multipliers 
(ADMM). ADMM is a special case of Douglas-Rachford 
splitting, [4], and it is related to other optimization algo- 
rithms, e.g. method of multipliers and Bregman iterative 
algorithms for £i problems, [5], [6]. For an overview of 
ADMM, we refer the reader to [7]. 

The most costly step in the proposed ADMM algorithm 
is the projection of an iterate to a set describing a feasible 
solution. We will show that this projection is equivalent to 
solving a Linear Quadratic (LQ) regulator problem with an 
additional linear term in the cost function. This problem can 
be solved efficiently using a Riccati recursion just as in [8]. 

We will apply ADMM to the recently introduced £i 
regularized Model Predictive Control (MPC), [9]. The £i 
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regularized MPC has an £{ regularized least-squares cost 
function. The motivation for £i regularized MPC is the re- 
duced actuator activity obtained when using i?i-norm penalty 
on changes of the input signal [9]. A detailed stability 
analysis of the closed loop system with £i regularized MPC, 
and results confirming sparse solutions are given in [9]. In 
£i regularized MPC, an optimization problem such as the 
one we consider is solved at each sampling instant. Hence, 
the sampling time puts an upper bound on the time that 
the optimization is allowed to take, and therefore efficient 
algorithms are needed. It is believed that ADMM is a 
preferred algorithm for this application based on the result for 
LASSO, [7]. We will see that this expectation is confirmed 
in numerical experiments. 

Also for traditional MPC an optimization problem has 
to be solved at each sampling instant, [10]. Because of 
this many different tailored optimization schemes have been 
developed to meet the real time requirements of MPC. 
Typically the optimization problem is a Quadratic Program 
(QP). There are mainly two different approaches that have 
been taken. One approach is to compute an explicit off- 
line solution to the QP which is stored in a look-up table, 
[11]. This facilitates very fast sampling, but is only feasible 
for small scale problems. The other approach is to com- 
pute the solution on-line as we propose, which is feasible 
also for medium- and large-scale problems. Among these 
approaches one can distinguish three different classes of 
methods: 1) Interior Point (IP) methods, [12], 2) Active- 
Set (AS) methods, [8] and 3) Fast Gradient (FG) methods, 
[13]. Riccati recursions play an important role also in IP 
and AS methods for MPC, [12], [8], since they can be used 
for these methods to efficiently factorize the matrix involved 
in the linear system of equations for the search directions. 
So far they have not been used for FG methods. For IP 
methods the Riccati recursion has to be re-computed for each 
iterate of the method. For AS methods it has to be updated, 
i.e. parts of the old solution can be reused but has to be 
modified. For ADMM it is possible to use the same Riccati 
recursion for all iterates. Computing the Riccati recursion, 
i.e. factorizing the matrix for the search directions, is the 
most time-consuming task for all these methods. However, 
the convergence performance is not the same for the different 
methods, i.e. it takes a different amount of iterations to 
reach a solution of satisfactory accuracy. For IP methods the 
number of iterations is typically 20-50 to reach very high 
accuracy. For active set methods the number of iterations are 
typically higher, however by considering gradient projection 
methods on the dual problem speed can be gained, [14], 
and similar results as for IP methods can be obtained. For 
fast gradient methods it has in [13] been shown how the 
number of iterates can be upper bounded to achieve a desired 



accuracy. Other recent relevant publications in relation to 
efficient methods for MPC include among others [15], [16], 
[17], [18] and the references therein. 

II. Control problem 

We consider an open-loop control problem of finding an 
input sequence that minimizes a finite-horizon cost function, 
given a model and an initial state. The problem is formulated 
as follows 
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where x, e R" is the state vector, «,- e R' is the input vector, 
yi E R'" and Zi E R'^ are auxiliary variables, and where 
||x||2^ = x^Ax. Formulation ([T]) captures the optimization 
problems that may occur in £i regularized MPC. For ex- 
ample, we can replace the input vector with the change of 
the input by augmenting the state vector and modifying the 
system matrices accordingly, see [10]. 

III. Alternating direction method of multipliers 
(ADMM) 

In this section, we provide a description of the key 
elements of ADMM. The description is a condensed version 
of the ones found in [19] and [7]. For a more rigorous 
overview, we refer the reader to [7]. 

A. Optimization problem 

ADMM is a numerical algorithm for solving optimization 
problems such as 



minimize f{x), 
subject to X E't^, 



(2) 



for some vector variable x G R", where f{x) is a convex 
function and '^ is a convex set. An equivalent problem to 
Ig) is 

minimize f{x)+I<^{xc), 
subject to x=Xc, 

where h-^{xc) is the indicator function of '^, [1]. 



(3) 



B. Augmented Lagrangian 

The augmented Lagrangian of optimization problem Q is 
defined as 

Lp{x,Xc,Xd)^fix)+I^f{x,) + {p/2)\\x-x,+Xd\\l, (4) 

where X/j is the dual variable corresponding to the equality 
constraint x = Xc scaled by 1/p, and p > is a tunable 
parameter There is no simple way of finding the optimal 
p, however there are guidelines in [7]. 



C. ADMM steps 

The ADMM algorithm consists of three main steps at each 
iteration k. The three steps are 

x'+' :=argmin{/(x) + (p/2)||x-x^+xf,||2} (5) 
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where n<^(x) denotes the Euclidean projection of a vector x 
onto a set '^. In the first step (|5]l, we minimize the augmented 
Lagrangian Q with respect to x, keeping Xc and x^i fixed. In 
the second step (|6]l, we minimize the augmented Lagrangian 
Q with respect to Xc, keeping x and x^/ fixed. In the third and 
last step Q, we update the scaled dual variable Xd- We then 
repeat all three steps until convergence. For more details and 
a complete convergence analysis, we refer the reader to [7]. 

D. Stopping criteria 

The ADMM algorithm is iterated until some stopping 
criteria are fulfilled. We use criteria based on the primal and 
dual residuals of the optimization problem. The primal and 
dual residuals of ([3| are 

e^, = ix'~4), 4 = -P(-^'-4"^)- 
We terminate the algorithm when 

Il4l|2<v^e*^ + e'-'='max{||x^||2,||xj||2}, 

where e*** > and e"^' > are absolute and relative toler- 
ances, respectively. For more details, see [7]. 

E. Over- relaxation 

We may use over-relaxation to improve convergence of 
the ADMM algorithm, [4], [20], [21]. When using over- 
relaxation we modify the update x^+' with 
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ax 
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where 1.5 < a < 1.8, in the second and the third ADMM 
steps, (|6]l and Q. For more details, see [7]. 

IV. Problem formulation and method 

In this section, we describe how the considered optimiza- 
tion problem in ([T]i can be solved using ADMM. 

A. ADMM formulation 

The optimization problem in ([T]i is on the same form as 
the optimization problem in (|2|. The vector variables are 

x=:(xo,...,x//), y^{yo,---,yH-i), 
u = {uq,...,uh-i), z = {zo,-..,zh-i), 

the objective function is 

H H 

/(x,y,M,z)=||x//|||g + ^||y,-_i||^+A^||z;_i||i, 

1=1 !=1 

and the constraint set is given by 

"^ ^ {(x,y,u,z)\xi =Ax,_i +Bm,_i, /=!,...,// 

Zi-i =Exi-\ +Fui-i, /= 1,...,//}. 



Thus, the AD MM formulation of optimization problem in 
Qis 



where 



minimize f{x,y,u,z) +I'^{xc,yc,Uc,Zc), 
subject to X ~Xc, y — yc, u = u^, z = Zc- 



(9) 



B. Step 1 ofADMM 

The first ADMM step, (|5]), is almost the same as the 
one for ii mean filtering in [19]. We solve 4i/+ 1 separate 
minimization problems because /(x,^, m,z) is separable in 
its arguments. For the vector variables x, y and u the 
minimization problems have a quadratic cost function and 
no constraints. The solutions are 
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x' / = 0,...,//-l, 



4+' = (2G+p/«r'p(4,//-4,//), 

yf+' = (2 + pr'p(3'L-^), '=0,...,i/-l, 
Mf+' =<,-<„ i = 0,...,H-l, 

where /„ denotes the identity matrix in R"^". For the vector 
variable z, the minimization problems are 



,<:+! 



:argmin{A||z,||i + (p/2)||z,~zJ,+4l|2}, (10) 



with component-wise solutions (zf^')j = -5^A/p((4.i~4 /)))' 
for i — 0,...,H —1, and j = l,...,p, where =^A/p denotes 
the soft thresholder operator, see [7]. 



C. Step 2 ofADMM 

The second step of ADMM, (|6]), consists of a projection 
of the vector 



(^k k k ^k\_(k+l.k .,k+l , k ,.k+l:,.k -k+l.k-, 

onto the constraint set '^, i.e. 

The projection can be formulated as the optimization prob- 
lem 

minimize |](xJ+\3;J+\mJ+i,zJ+1) - (x^,3;^,4,z^)|l2, 

subject to (4+i,3;J+i,M^+Sz*+i) G '^. 

(11) 
To simplify notation in the rest of this section, we will drop 
the use of super script k and k+l. An equivalent optimization 
problem to the one in ( [TT| is 



minimize v^ ^v + q^ v 
subject to ^v = g, 



(12) 



V — {Xcfl,Uc,0,---,Xc,H-l,Uc^H-UXc,H), 
g^ (Xc,o,0,...,0), 

q= {ro,SQ, ...,rH-uSH-uO), 
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The symbol (E) denotes the Kronecker product. The optimiza- 
tion problem in ( [T2] i is an equality constrained minimization 
problem. As such, its solution is equivalent to the solution 
of its Kamsh-Kuhn-Tucker (KKT) conditions, [1]. The KKT 
conditions of the optimization problem in ( [T2j i are 



2^ 







-^ 



= 0, 



(13) 



with w — (v, v^/), where \\i is the Lagrange multiplier cor- 
responding to the equality constraint ^v = g. The KKT 
conditions in ( [T3] l are a system of linear equations and can be 
efficiently solved using a Riccati recursion as described in the 
Appendix. The solution to the optimization problem in ( [TT| ) 
is obtained by extracting Xc and Uc from v, and calculating 
yc and Zc from the equations defining the constraint set "^ . 

D. Step 3 ofADMM 

In the third ADMM step in (|7]), we update the scaled dual 
variables, i.e. 

V -"dj + y'^i -«6-,; )' 

^d,i ~ ^d,i ^ \^i ^c,i I ' 

V. Example 
In this section, we describe how the ADMM algorithm. 
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proposed in IV can be used to solve an £i regularized MPC 
problem without inequality constraints. 

A. Model 

We consider a linear and discrete model of the plant. The 
model is given by 



x{t + \) =Ax{t)- 
y{t)=Cx{t), 



-Bu{t), 



(14) 



where x{t) G R" is the state vector, u{t) e R' is the input 
vector and y{t) E R™ is the output vector. 



B. Cost function 

The control objective is to drive the output vector to zero, 
namely the regulator problem [10], while using a piece- 
wise constant input signal. Such a control objective can be 
described by the cost function 



y(0 = p(f+//,|0l|2g+i£\||y(f + /-i|0llie- 



^lf='ll|AM(f + !-l|0 



(15) 



The cost function penalizes the terminal state, output de- 
viation from zero and non-constant input signals. In ( [T5] l, 
x{t + i\t) and y{t + i\t) are the predicted state and output 
vectors, respectively, at time t + / given measurements up to 
time t and the model in ( [T4j l. Moreover, 

AM(f + /|f) = u{t + i\t) - u{t + i - l|f), 

where u{t-\-i\t) is the predicted input vector given measure- 
ments up to time t and the model in ( [T4| i. The prediction 
and control horizons are denoted Hp and //„ respectively, and 
we assume that ls.u{t + i\t) = for all / > //„. The matrices 
2, 2 G -^nxm ^jjj j-jjg scalar A are weights. We require that Q 
and Q are positive semidefinite, and that A is non-negative. 

C. Optimization problem 

The control objective can be achieved by minimizing the 
cost function in ( [T5| ) given the model in ( [T4] i in each time 
step t, in accordance with the receding horizon idea, [10]. 
We can formulate the optimization problem as 

minimize V{t)^ 

subject to x{t + i\t) ^Ax{t-\-i- l\t) + Bu{t + i - \\t), 
/= l,...,Hp, 
x{t\t)^x{t), 

y{t + i-l\t)=Cx{t + i-l\t), i=l,...,Hp. 

(16) 
The optimization problem in ([T6| is similar to standard for- 
mulations as the one found in [10]. The significant difference 
is the use of the £i-norm of Aii{t + i\t) instead of the £2- 
norm in the cost function in ( fTSJ ). The former typically 
promotes sparse Au{t + i\t) for / = 0, . . . ,//„ — 1, while the 
latter promotes small but non-zero elements of Aii{t + i\t) for 
i — 0,... ,Hi, — 1, [1]. To simplify notation in the rest of the 
paper, we denote x{t + i\t) as xi, y{t + i\t) as yt, and so forth. 



E. MFC formulation 

We consider the optimization problem in ([16]). We set the 
predicted output vector to be the predicted state vector, and 
the prediction horizon equal to the control horizon, that is. 



minimize ||x//||2 = +52fLi \\xi-i\\2,Q + ^LiLi I|Am, 



subject to Xi — Axi^ 1 + Bm,_ 1 , i — l,...,H. 

(17) 
The optimization problem in ^Tf) can be reformulated by 
replacing m, with Am,- in a similar way as in [10]. We 
introduce three new vector variables i, e R"+', y, G R" and 
Zi G R' in the following way 

minimize H^Hj g + lfLi llyz-illi + '^L^Li lk;-i|li, 
subject to Xi ~ Axi^ 1 + BAm,_ 1 , i = l,...,H, 

yi-i = Ci,_i +DAm;_i, /=!,... ,i/, 
Zi-i =Exi^\ +FAm,_i, / = 1,...,//. 

(18) 
Here, x; is the state vector augmented with the input vector 
at the previous time step, i.e. x; — {xi,Ui^i). The matrices A, 
B and C are given by 





'A B ' 




■ B ' 


A = 


. // . 


, B = 


_ h _ 



, c=[ 







where c is chosen such that Q = c^c, and the matrices D, E, 
F and Q are given by 



£) = 0, £ = 0, 



Q 



Q 




The optimization problem in ( [T8[ l is of the same form as the 
one in ([TJ, with Am, acting as the input. Therefore, we can 
solve ( fTSJ ) efficiently using ADMM and a Riccati recursion as 
described in Section IV Note that we in this particular case 



can pre-calculate Fi, Hi, Gi, and 5, in the Riccati recursion, 
before we start the MPC iterations. 

VI. Numerical examples 

In this section, we apply the ADMM algorithm on an £1 
regularized MPC problem. All the examples are performed 
with p = 1. To improve convergence we use over-relaxation 
with a = 1.8 and we warm-start each ADMM iteration with 
the variable values obtained in the previous MPC iteration. 
We use stopping criteria ([8| with e"''* = 10"^ and e'''"' = 
10-4. 



D. Receding horizon 

The optimization problem in ( [T6] l is solved with respect to 
the input vector m,- for / = 0, ...,//„ — 1. The input vector at 
the first time step, mq, is applied to the plant. The state vector 
is updated according to measurements and, if necessary, an 
observer. The optimization problem in ([T6| is updated and 
solved again. The described procedure is repeated until some 
final time step. Note that closed loop stability cannot be 
guaranteed for all values of A, see [9]. Typically, a terminal 
cost penalty is used to obtain closed loop stability, if possible. 



A. Example 1: Quadruple water tank process 

I) Plant: The plant is the quadruple water tank process 
presented in [22]. The process is shown in Figure [T] where 
x= (xi,X2,X3,X4) are the water levels, u = (mi,M2) are the 
pump voltages, and 71 = 72 = 0.625 are the parameters asso- 
ciated with the valves. The area of the cross-sections of the 
outlets of each tank are (01,02,^37^4) =(0.17,0.15,0.11,0.08) 
cm^, the area of the cross-sections of each tank are A = 15.5 
cm^ and the parameters associated with the pumps are ki = 
fe2=4.14cm3/(sV). 




Fig. 1. Quadruple water tank process. 



2) Model: We obtain a linear model of the process by 
linearizing the nonlinear plant description given in [22] 
around its equilibrium points. The linearized model is 
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The equi- 



librium points of the plant are x^ = (15,15,3,12) cm and 
u^ = (7.8,5.25) V. The linear model is discretized assuming 
zero-order hold sampling at a sampling rate of 1 Hz. 

3) Simulation: We set // = 5, Q = h, and 2 = 0. The plant 
is initialized with ji:(0) = (16, 16,4, 13) cm and m(0) = u^. A 
Kalman filter is used to estimate the complete state vector 
during simulation. The MPC iterates for 10 time steps. We 
perform the same MPC simulation for X equal to 0.05, 0.1, 
2 and 5. The applied input sequences are shown in Figure 
[2] and the output sequences are shown in Figure [3] We see 
that the applied input signal varies over time for low values 
of A. As A gets larger, the input signal becomes piece-wise 
constant, and eventually completely constant. We also see 
that a more restrictive control strategy, i.e. a high value of 
X, gives worse control performance in terms of response time 
and static error. 

B. Example 2: Number of iterations in ADMM 

Figure |4] shows the number of iterations required in 
ADMM for fulfilling the stopping criteria in Example 1 . We 
also investigated the number of iterations required without 
warm-starting the algorithm. The conclusion is that a warm- 
start improves the convergence of ADMM when the plant 
inputs are close to constant and no rapid changes in the 
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Fig. 2. Appl ie d inpu t seq uences. The applied input sequence denoted 

\ ) , \ ) , {■-■-) and j I conesponds to X equal to 0.05, 0.1, 2 

and 5 respectively. 
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Fig. 3. Measur e d outpu t seq uences. The measured output sequence denoted 

\ ) , \ ) , \- - ■ -) and \ \ c orresponds to X equal to 0.05, 0.1, 2 

and 5 respectively. The sequence )°|°|°) corresponds to the equilibrium point 
of the water levels. 



plant states occur. When this is not the case, we get similar 



performance with and without warm-start. It is natural that 
the benefit of warm-starting is greater the less the states 
move. 




Fig. 4. Number of iterations in ADM M requi re d for f ul filling the st opping 

criteria. The iteration sequence denoted \ \ , { \ , {- -) and j ) 

corresponds to A equal to 0.05, 0.1, 2 and 5 respectively. I'he number of 
iterations required drops when no rapid changes in plant inputs or states 



C. Example 3: Convergence of ADMM 

Here we investigate the same set-up as in Example 1. We 
only consider A =0.1 and the first optimization problem 
solved in the MPC iterations. We calculate the error of the 
cost function in ([TSj for each iteration in ADMM. The error 
is defined as e'^ ^ V* — V'', where V* is the true optimal 
value of the cost and V^ is the value obtained in ADMM 
iteration k. The true optimal value is approximated with the 
solution obtained by running ADMM for 1000 iterations. The 
optimal value is verified using CVX, a package for specifying 
and solving convex optimization problems, [23]. CVX calls 
the generic SDP solvers SeDuMi [24] or SDPT3 [25] to 
solve the problem. We choose to use SDPT3. The resulting 
error is shown in Figure |5] The true optimal value is V* = 
4.58108, and the final value obtained from ADMM is V^^ = 
4.58105, where 264 is the number of iterations required to 
fulfill the stopping criteria. A rapid drop in the error occur in 
the first iterations in ADMM. The ADMM algorithm iterates 
until the stopping criteria are fulfilled, however, for improved 
visibility of the drop we only show the first 50 iterates. Note 
that since the ADMM solution is not necessarily feasible it is 
possible to achieve a value of the cost function at iteration k 
that is lower than the optimal one. The corresponding primal 
and dual residuals are shown in Figure [6] We see a rapid drop 
in error and residuals for the first 20 iterations in ADMM 



(e-"" = -0.03, e;" = 0.11 and ef = 0.06), confirming that 
ADMM converges fast to a moderate accuracy. 

D. Example 4: Time of iterations in ADMM 

We consider the set-up in Example 1 with A = 0.1 and a 
prediction horizon H varying from 5 to 100 in steps of 5. 
We only consider the first optimization problem solved in 
the MPC iterations and we fix the iterations in ADMM to 
1000. We calculate the mean value of the time required for 
an iteration in ADMM. Figure 17] shows the resulting means 
with respect to the prediction horizon. We see that the mean 




Fig. 5. En'or of cost function. The eiTor for each iteration k in ADMM is 
shown. 



-a 
e 



10 



10 



ion 


1 


^^--^=-^:^ 




0-5 - 


^^v 


-10 - 


\v 


-IS 


1 1 1 ^)>w 



200 



400 600 
k 



800 1000 



ration of 



Fig. 6. Primal and dual residuals. The primal \ ) and dual 

residuals are calculated for each iteration in ADMM in the first iteration ( 
MPC. The vertical line shows where the stopping criteria are fulfilled. A 
rapid drop in the residuals occur in the first 20 iterations in ADMM. 



time of the iterations in ADMM is linear in the prediction 
horizon. This is expected since the computational cost of the 
Riccati recursion is linear in H, [26]. 




Fig. 7. Mean value of time required for an iteration in ADMM. The 
mean values are calculated over 1000 iterations in ADMM, for a prediction 
horizon H varying from 5 to 100 in steps of 5. The mean time is linear in 
the prediction horizon. 



E. Example 5: Required accuracy 

Example 3 shows how close the ADMM solution is to 
the optimal one for e"'" — 10^^ and e""' = lO^'* in stopping 



criteria (|8]l. However, the stopping criteria used may be to 
conservative with respect to required control performance. 
For example, if we in Example 1 with A = 0.1 restrict 
the number of iterations in ADMM to 10, we can have a 
sampling rate of 100 Hz in the MPC, see Figure IT] The input 
signals obtained with both 10 and 1000 iterations in ADMM 
are shown in Figure [8] The corresponding output signals are 
shown in Figure |9] We see that, although the signals differ 
from each other, they still have the same over-all behavior. 
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_ Ap plied i nput sequences. The applied input sequence denoted 

j- - - ( and \ ) con-esponds to 10 and 1000 iterations in ADMM 

respectively. The sequences differ, however the over-all behavior is the same. 



VII. Conclusion 

We have derived a method for solving optimization prob- 
lems with an l\ regularized cost function subject to recursive 
equality constraints. The optimization problem occurs in 
control applications, e.g. i\ regularized MPC. The method 
is based on the ADMM algorithm. We have showed that 
the costly projection step in ADMM is equivalent to solving 
an LQ regulator problem with an additional linear term in 
the cost function. Such problems can be efficiently solved 
using Riccati recursion. Future work consists of expanding 
the proposed method to (,\ regularized cost functions subject 
to both recursive equality and inequality constraints. 

Appendix 
Riccati recursion 

We use a Riccati recursion to solve the projection problem 
(|6]l, as in [8]. We showed in Section [IV-C| that the solution 
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Fig. 9. Meas ured ou tput sequences. The measured output sequence denoted 
and { ) con'esponds to 10 and 1000 iterations in ADMM 



respectively. The sequences differ, however the over-all behavior is the same. 

to (|6]l is equivalent to the solution of a system of linear 
equations. 
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where =S and s^ are block-diagonal matrices defined as 

Qn Qn 



^- 



Q 
Q 



Q 



and 



^ = 



7 
-A -B I 



Oil Qii 



-A -B I 

The vectors ^, A, rs and rx can be divided into sub-vectors 

t, = (xo, MO, •••,«//-!, -«//), A = (Ao,...,A//), 

ri, = ('•x,0,'-»,0,---,'-»,//-l,rT,//), a = ('"A,0,---,'"A.//), 

where x,- is given by the system equations 

x;+ 1 = Axi + Bui + rx .,■+ 1 . 

It is shown in [26], that there exists a matrix 5, and a vector 
•P,- such that 

A,+5,x,=^,, i^O...H, 

where Sh — Q, ^h = f.x.H, and xq = r^Q. The matrices 5, and 
vectors ^; for i — Q,...,H can be found through backward 



recursion. We then obtain ^ and A through forward recursion. 
The algorithm is as follows [26]: 
Backward recursion: Update Si and 'P,, 

^■+1 = Qn +A^5,+iA, Hi+i = Qn +A^S,+iB, 
Gi+i = Q22+B'^Si+iB, \^i+i=^>i+i-Si+irx^i+i, 

"Vi = r^s+A^Wi+\ -Hi+iG;_^\{r^j + B^\i/i+i). 
Forward recursion: Update A;, m,+i and x,+i, 

A. Unstable model 

If A is unstable, the Riccati recursion might not provide 
the correct solution to ( [T9] l. To avoid this, we pre-stabilize 
the state-space equations using state feedback control, see 
[26], [27], [28]. That is, we let 



(20) 



x,+ i = Ax,- +Bui, 
Ui = -Lxi + Vk, 

where L is the feedback vector. We can reformulate ( |20] i as 

Xi+i = {A-BL)xi+Bvi, 

and treat v,- as the unknown input signal. The solution 
obtained from the Riccati recursion will be the values of 
Xi and V,. The solution in terms of the original parameters x, 
and Ui, can be obtained as 
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